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Distance estimation with efference copies and optical 
flow maneuvers: a stability-based strategy. 
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Abstract 

The visual cue of optical flow plays a major role in the navigation of flying insects, 
and is increasingly studied for use by small flying robots as well. A major problem is 
that successful optical flow control seems to require distance estimates, while optical 
flow is known to provide only the ratio of velocity to distance. In this article, a novel, 
stability-based strategy is proposed to estimate distances with monocular optical flow 
and knowledge of the control inputs (efference copies). It is shown analytically that 
given a fixed control gain, the stability of a constant divergence control loop only de¬ 
pends on the distance to the approached surface. At close distances, the control loop 
first starts to exhibit self-induced oscillations, eventually leading to instability. The 
proposed stability-based strategy for estimating distances has two major attractive 
characteristics. First, self-induced oscillations are easy for the robot to detect and 
are hardly influenced by wind. Second, the distance can be estimated during a zero 
divergence maneuver, i.e., around hover. The stability-based strategy is implemented 
and tested both in simulation and with a Parrot AR drone 2.0. It is shown that it 
can be used to: (1) trigger a final approach response during a constant divergence 
landing with fixed gain, (2) estimate the distance in hover, and (3) estimate distances 
during an entire landing if the robot uses adaptive gain control to continuously stay 
on the ‘edge of oscillation’. 


1 Introduction 

A major challenge in robotics is to achieve autonomous operation of tiny flying robots 
such as 25g “pocket drones” m or more extremely the 80mg Robobee [31] . Flying insects 
provide a rich source of inspiration for solving this challenge, since they are able to navigate 
successfully on the basis of a very limited sensory and processing apparatus. Flying insects 
rely heavily on optical flow, i.e., the apparent motion of world points caused by the relative 
motion between the observer and the environment m Also for robots this cue is promising 
as it can be extracted from a single passive camera, implying light weight and relatively 
little energy consumption [Huang Eg. 
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It is well-known that optical flow can convey information on the ratio between distance 
and velocity, but that additional information is necessary to disentangle these quantities. 
This ‘scaling’ can be performed by means of maneuvers that have a known range or velocity 
(as in parallax m. or with any sensor that directly or indirectly provides distance or 
velocity estimates, e.g., accelerometers or sonar mm- 

Flying insects such as fruit flies do not seem to have scaling sensors [IS]. In addition, it 
is unlikely that they know the range of their maneuvers since they only have access to their 
air speed, not their ground speed |B]. Nonetheless, flying insects are able to navigate suc¬ 
cessfully. They follow straightforward strategies based on optical flow observables. First, 
it was found that honeybees keep ventral flow (defined as d x = ^f ) constant when they 
perform grazing landings [39]. This strategy was studied for flying robots and simulated 
spacecraft [311 SSj 121 E2S] . However, these studies required a “pitching” law, ultimately 
depending on an estimate of the height and velocity. The main reason seemed the absence 
of control of the vertical dynamics. Therefore, more recent landing strategies also include 
the flow divergence ( D , typically defined as with the positive z, v z axis pointing up) 
or time-to-contact (r = F 2 ) [22| [2T, EE Ij- These strategies either enforce a decreasing 
time-to-contact (as humans do for car braking [2B]) or keep the divergence constant (as 
honeybees do for straight landings 0)- Recently, in [T5], ventral flow and divergence were 
combined with slope estimation and an anemometer in order to fly over uneven terrain, 
without a need for inertial frames or scaling sensors. 

The absence of velocity and distance information has been a defining property of the 
work on optical flow based navigation strategies. However, some observations suggest that 
for successful navigation distance estimates are still of importance. First, for successful 
optical flow landing control, the gains of the controller are tuned to a range of specific 
heights and velocities. For instance, for the constantly decreasing time-to-contact landings 
in m , the gains are scheduled by means of r, along a trajectory starting at a specific initial 
height and velocity. Second, while for a constantly decreasing time-to-contact landing it is 
possible to perform a final landing procedure (when the time-to-contact is almost zero), it 
is not obvious how such a procedure should be triggered in a constant time-to-contact or 
divergence landing. 

It is actually possible to estimate distances based on optical flow, if a robot has access 
to copies of its control input signals (referred to as “efference copies” in the biological 
literature). Two main strategies have been proposed so far in the literature. 

The first strategy for estimating distances is well-established in the field of Image Based 
Visual Servoing (IBVS). In the task of IBVS a robot has to regulate certain target features 
to desired image locations, without knowing a priori their dimensions or geometry. Quite 
early it was observed that if a fixed gain is used to map an image error to a control input, 
the response of the system depends on the distance to the target [3]. Given a model of the 
robot’s dynamics, the effects of the robot’s actions on the image features can be used to 
estimate the distance to the target in an adaptive gain control scheme (see, e.g., 0). 

The second strategy for estimating distances was recently proposed by [6J. They showed 
that during a constant divergence approach of an object, the control input u (the thrust) is 
proportional to the distance d - and that efference copies can then be used as a stand-in for 
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distance. This solution was tested on a camera mounted on a rail and successful distance 
estimates were obtained by integrating information during the entire approach. 

Both strategies above assume that the robot’s accelerations are completely determined 
by the control inputs, ffowever, flying robots and insects are also subject to significant 
external accelerations caused by wind or wind gusts. 

The main contribution of this article is the proposition of a novel strategy for distance 
estimation with optical flow and efference copies. The central idea is to estimate distance 
by exploiting the self-induced oscillations that result from the fundamental imperfection 
of fixed-gain optical flow based control. The strategy provides two major advantages with 
respect to the previously proposed strategies: (1) it works straightforward even in the 
presence of wind, and (2) it can be used in a broader range of maneuvers, and notably in 
hover. 

The remainder of the article is organized as follows. In Section^ the distance estimation 
method previously proposed in [6j is explained and investigated in the presence of wind. 
Subsequently, the new stability-based strategy to distance estimation is proposed in Section 
[3] In Sectioned is shown how self-induced oscillations can be used to trigger a final landing 
procedure. Subsequently, in Section 0 adaptive gain control is introduced to determine the 
distance in hover and to estimate distances along an entire constant divergence landing. In 
Section 1 the hirelings are empirically verified with real-world experiments using a Parrot 
AR drone. Finally, the results are discussed in Section [7| and conclusions are drawn in 
Section [HJ 

2 Thrust-based Distance Estimation 

In this section, the distance estimation method proposed in [6] is explained and investigated 
in the context of a (robot) landing task. For reasons that will become clear below, this 
strategy will be referred to as a ‘thrust-based’ method. It is shown that the application of 
this method is complicated by the presence of wind. 

Figure Q] shows the definition of the axes as used in the formulas. In order to keep the 
equations as uncluttered as possible, most of the analysis in this article focuses purely on 
the z-axis. Generalizations to other looking and movement directions are easily made (see, 
e.g., Appendix O, as are the inclusion of different attitudes or a displacement and rotation 
of the camera with respect to the body center of mass. 

Importantly, as ‘visual observable’, the variable d z is introduced, which is the relative 
velocity —. It is related to the divergence as = —D, and hence a constant or zero diver¬ 
gence landing is also a constant or zero tb landing. The relative velocity d 2 is determined 
on the basis of the spatial gradient of optical flow in the image, but the algorithms to do 
this are considered outside the scope of the article and the reader is referred to, e.g., ED- 
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Figure 1: Axis system employed in the article. 

2.1 Relation between u and z 

As was shown in [6], when perfectly following a constant divergence landing, there is a 
straightforward relation between the control input u z (thrust) and z. Let us start from the 
visual observable '() z : 

** = - ( 1 ) 


If ti z is differentiated over time, it results in: 


Using Eq. [D 


A = ^ 

z z 2 


z = 


til + ti z 


( 2 ) 


(3) 


So, z is expressed in terms of ‘observables’ (ti z and its time derivative ti z ), and the vertical 
acceleration a z . Adding an accelerometer to the robot will allow for estimation of z (e.g., 
[23]). However, it requires using the time derivative of ti z , which typically induces a lot of 
noise in the observation. 

Van Breugel et al. [6j describe the equations of motion in a linear state space model 
as: 


i(f) 

v z (t) 


0 1 
0 0 


z(t) 

V z (t ) 


u z (t) 


This model entails that u z = a z . Making this assumption, Eq. [3] becomes: 

Vj Z 


z = 


til + ti z 


(4) 


(5) 
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In [6] it has been investigated what would happen if a perfectly constant divergence landing 
is performed, $2 = 0 and d z = —c 2 (c 2 to indicate that d z has a negative set point): 


Eq. IH] shows that during a perfect constant divergence landing, the thrust u z is a scaled 
version of the height. Hence, it can be used as a stand-in for the height. For this reason, 
the method is referred to in this article as a ‘thrust-based’ method. 


2.2 Dealing with wind 


The assumption a z = u z effectively implies that any possible gravity is subsumed under 
the w 2 -term: 

v! 

a z = ~g H —~ = u z (7) 

m 

, where u' z = m(a z + g ) is the actual upward thrust generated by the robot in Newton. 
Two observations are in place. First, it is common in control system theory to subsume 
factors such as gravity and the mass into u z , as the calculation of u' z is straightforward 
and does not depend on any state variables. Second, it is unlikely that a robot’s command 
signal (and hence efference copy) is equal to u' z . A command signal u" (in any unit, e.g., 
the commanded RPM of a robots’ propellors or the flapping frequency of an insect) will 
typically have to be mapped to u' z (in Newton) with an actuator effectiveness estimate. 
If such a mapping is unknown, the estimated height z from Eq. [H] will have a possibly 
nonlinear relation to z\ 

z = g(z). (8) 

As long as the function g is invertible, this mismatch in unit essentially does not matter. 
Taking this reasoning into account, for the remainder of the article a copy of u' z will be 
referred to as an ‘efference copy’. 

Eq. [7] is only valid if the control force is the only force acting on the robot. In a vacuum 
environment this can be approximately correct, e.g., for a moon landing. However, a robot 
flying in the air will undergo additional accelerations. Most important is the drag force, 
which depends on the robot’s movement relative to the air surrounding it: 


f D = sign (u air ) - pC D Av, 


2 

air’ 


(9) 


where: 


^air ^wind 


( 10 ) 


and sign(u air ) indicates the directionality of fo along the ,3-axis. This leads to an additional 
acceleration: 


a z = ~g + 


< + Id 


m 


(ii) 


, where fo is a time-varying value involving a non-linear function of v z and an uncontrolled 
variable u win a. 
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It is informative to study the effect of drag and wind on the required u' z , when the 
robot follows a perfect constant divergence landing. In such a landing, = —c 2 at every 
point by definition. The acceleration resulting from drag and the required u’ z are calculated 
to match. Figure [2] shows the value of u' z for a robot with assumed values of m — 1kg, 
p = 1.204, and g = 9.81m/s 2 . Furthermore, Cd = 0.25, and A = 0.25, which are rather 
conservative values representing for instance a small quad rotor. The figure shows the 
thrust u' z versus height z in different environmental conditions: vacuum (red solid line), 
wind still (black dotted), in a downward wind of — lm/s (blue dashed), and in an upward 
wind of lm/s (green, dashed-dotted). 



Figure 2: Thrust u' z versus height z in different environmental conditions (see legend). Two 
black lines illustrate an estimation error when the wind is unknown. The lines indicate 
that the thrust with a downward wind of —lm/s at z — Orn corresponds to a thrust in a 
wind-still environment at a height of z ~ 4.5m . 

Three main observations can be made from Figure |2j (1) the red solid line in vacuum 
indeed indicates the derived linear relationship between u z and z, (2) in a wind-still envi¬ 
ronment, this relationship is nonlinear, but invertible - and hence as useful, (3) adding a 
modest wind speed to the equation already makes significant differences to u! z . Two black 
lines indicate that the thrust with a downward wind of —1 at z — 0m corresponds to a 
thrust in a wind-still environment at a height of z ~ 4.5m. The curve for a wind of lm/s 
does not even match that in a wind-still environment. Hence, even with a perfect constant 
divergence landing these differences distort the relationship between u z and z. 

So, in order to retrieve the right u z -curve (see Figure [2]), the wind will have to be 
measured. Flying insects may be able to measure the wind velocity with their hairs [IS] , 
while robots can be equipped with an anemometer, as in be Although measuring the wind 
velocity is a possibility, this does not lead to an easy solution of the distance estimation 
problem. First and foremost, controlling a constant divergence landing in a real system 
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means dealing with sensing and actuation delays, visual inaccuracies, and changing external 
factors such as wind gusts. Any real system will thus deviate from the perfect landing 
profile, leading to the necessity of command signal variation and hence a varying u' z . Indeed, 
this effect can be seen in the raw 5 estimates in |6j for the camera-system mounted on a rails. 
In [6] they tackled this issue with a robust estimation scheme that integrates information 
during the entire approach. The problems are much worse though for a freely flying system, 
which is also subject to wind. 


3 Instability of Constant Divergence Landings 

The onset of the work in this article came from the difficulty of making optical flow based 
landing strategies work on real Micro Air Vehicles (MAVs). The study in [26j is especially 
instructive, since multiple control laws are studied, of which the gains depend to different 
extents on the initial conditions. Especially important are the remarks and observations 
on the difficulty of getting the optical flow control to work close to the landing surface, 
even suggesting to disengage the optical flow control under a threshold distance (as is also 
done in, e.g., [42]). 

Please note that it is not easy to identify the cause of oscillations close to the surface 
when performing vision-based control with a real MAV. Oscillations close to the ground 
can potentially be explained by the “ground effect” due to the MAV’s downwash and by 
the vision process (more blurry images, larger optical flow vectors that are harder to track). 
Below it will become clear though that even without such additional effects, the instability 
problem arises. 

Based on the dependence of the gains on the actual height and velocity, in this article a 
novel strategy is proposed for distance estimation: stability-based distance estimation. In 
this section, first the fundamental reason for the gain tuning problem of optical flow based 
control is explained (Subsection 13.11) . Subsequently, a specific control law is assumed and 
it is shown analytically that the stability directly relates the control gain to the height 
(Subsection 13.2[) . 


3.1 Fundamental reason for gain tuning 


In order to see why a single gain is only stable for a given range of heights, let us go back 
to Eq. |3l and reorganize it as: 


ti z 



z 


( 12 ) 


Assuming a z as in Eq. ITU 


<+/d 


ti z = 


-til 


Differentiating this formula with respect to u' z gives: 


(13) 


dti z 1 
du' z mz 


( 14 ) 
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Eq. IT-11 shows that a change in thrust has a larger effect on the change in $ z close to the 
landing surface than far away (independently of external accelerations). The constancy of 
this relation depends on the constancy of the robot’s mass, which applies for landing quad 
rotors but not for spacecraft. A spacecraft will loose more and more mass when landing, 
hence aggravating the effect in Eq. fl4l 

Thus, a control gain that leads to a satisfactory control performance far away from the 
landing surface, will lead to very large effects on d 2 and hence indirectly on close to the 
landing surface. Although theoretically one could use directly for distance estimation 

(z = j n practice this value is extremely noisy, since it involves a partial derivative 

of the already noisy value i9 Z . 


3.2 Relation between control stability, gain, and height 

In this subsection, it will be shown that the larger influence of u z on d 2 at lower heights 
eventually leads to instability at a specific height. For the analysis, the following constant 
divergence control law is studied, as was used in |6j: 


u z = K z (&* e -0), 

where it is easy to see that in a noiseless, delayless system: 


Vj Z 


lirn — = (tf* - i) z ) = 0, 
I< z ^oo A 


(15) 


(16) 


and hence i? z = d* z . Real systems always have some noise and delay. Below, it will be shown 
that just discretizing the control with Zero-Order-Hold (ZOH) will already introduce an 
instability into the system. 

Let us assume the state space model of Eq. [H ignoring drag for the moment. The 
corresponding observation is: 

y = *. = -, 


z 

which is a nonlinear function. Linearizing the state space model gives: 

Ay(t) = ±][Az(t) Aia(f)] T , 

so that the state space model matrices are: 


A = 


0 1 
0 0 


B = 


C=[^ i] O = (0], 


(17) 


(18) 


(19) 


As mentioned, the continuous system without noise and delay is stable (Eq. fTBl) . 
Here, the discretized system is studied, which has the following state space model matrices 
corresponding to the continuous ones in Eq. [19j 


$ = 


1 T 
0 1 


r = 


rp 2 ' 
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Figure 3: Root locus plots of a ZOH-model with T = 0.03s for z = lrn and z = 10m. 
The plots are equal in shape, but for each point K z is different. This is indicated for the 
(approximations to) the values of K z at which the system becomes unstable (w = —1). 


where T is the discrete time step in seconds. The transfer function of the open loop system 
can be determined to be: 

G(w) = C(wl - $) -1 r, (21) 

{zT - \v z T 2 )w -zT - \v z T 2 
z 2 {w — l) 2 

where we use w as the Z-transform variable, since z already represents height. The feedback 
transfer function is: 



G{w) 


K z ((zT - \v z T 2 )w -zT- \v z T 2 ) 
z 2 (w - l) 2 + K z ((zT - \v z T 2 )w -zT - \v z T 2 )' 


(23) 


Eq. [22] shows that, given a gain K z and time step T, the dynamics and stability of the 
system depend both on the height z and velocity v z . 

Figure [3] is the root locus plot of G(w). For the Z-transform any pole in the unit circle 
is a stable pole. The different elements of the plot can be found back in Eq. [23] The 
numerator indicates a single finite zero at: 


w 0 = 


zT + \v z T 2 
zT — \v z T 2 ' 


(24) 


Given z > 0, T > 0, and v z < 0, wq is positive and due to the typically small T values 
slightly smaller than 1. There is also a negative infinite zero. Since the denominator of 
Eq. [23] is a second order equation in w, G{w) has two poles, both of which are stable. For 
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K z = 0, the poles are located at w = 1. As the gain K z increases, one pole moves toward 
the finite zero and one toward the infinitely negative zero. This latter pole is interesting, 
because as soon as it passes w = —1, the system becomes unstable. Setting w = — 1 in the 
denominator and equating it with 0 gives an equation for the unstable gain value K z : 

Kz = (25) 

This implies that given a specific gain K z and time step T, there is always a height at 
which the control system will become unstable: 

z = \k z T. (26) 

When using a fixed gain, a tradeoff has to be made between the control performance at a 
larger height (requiring a large gain) and at lower heights (requiring a small gain). Please 
remark that the instability at lower altitudes is exactly what is observed for robotic systems 
performing constant divergence landings! 


3.3 Including wind 


In this subsection it is shown that also for a model with wind, the instability of the system 
given a K z depends on z. Including wind leads to a changed formula for v z (see Eqs. [9] 
and |TT]) : 


v z =u z + sign(u wind - v z ) 


1 

2m 


v z ) 


2 


(27) 


, where from now on (5 will stand in for the constant ApC D A. In order to obtain a linear 
state space model, this equation is linearized to obtain: 


Aib = A u z - (sign(u wind - v z )fi(v wind - v z )) Av z (28) 


, where the factor multiplied with Av z is a constant in the linearized model (u w i nd and v z 
are given the values at the linearization point). In order to avoid cluttering the formulas, 
this constant is represented by p = sign(u wind — v z )/3(v win d — v z ) (where p > 0). This leads 
to the following continuous, linear state space model: 


~P 


B = 


A = 

which results in the Z-transform matrices 

$ = 


c=h? a o=io], 


r (l-e- pT )l 

P nn 

r - 

0 e~ pT 



T _ (1—e~ pT ) 
P P 2 

(l-e~ pT ) 
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c=i^ a b=io], 


(29) 


(30) 

(31) 


10 













and a rather involved transfer function H(w). The root locus plot of H(w ) is very similar 
to that of G(w), also with a negative infinite zero. Following the same procedure as for 
the ‘vacuum’ model, setting w = —1 and equating the denominator with 0, gives: 

= _ (2 P 2 + 2 p^T)^ _ 

( 2e pT — 2 — Tp — Tpe pT )v z + (2 pe pT — 2 p)z' 

Eq. [33] includes many instances of p that depends on v z and n w ; n d, and also contains a term 
in the denominator with v z . This suggests that there is no fixed linear relation between 
K z and z. However, closer inspection shows that the term (2e pT — 2 — Tp — Tpe pT ) ~ 0. 
Rearranging terms, this leads to: 


2 p 2 + 2 p 2 e pT 
2 pe pT — 2 p 


(33) 


which still depends on p. It turns out that the fraction that is multiplied with z is almost 
identical to the fraction Indeed, solving Eq. [33] for different variable settings gives 
almost identical results to Eq. [23] 


3.4 Continuous system with drag and delay 

The previous subsections discussed discretized models for which unstable K z exist, and 
where the limit case can be analytically expressed as a function of z. In order to show 
that continuous systems encounter the same type of phenomenon, Figure® shows the root 
locus plot of a continuous system with wind, a delay of At = 0.03s and c 2 = 0.01. The 
delay introduces two zeros in the (unstable) right-hand plane of this continuous root locus 
plot. A few values of K z are plotted, where the poles get an imaginary component and 
where they cross the imaginary axis (shown for z = 10m and z = lm). 


4 Stability-based Distance Estimation 

In the previous section, it was established that the instability of a constant divergence 
landing system depends on the height. If this instability can be detected by a robot or 
insect, it can be used for distance estimation! 

4.1 Detection of Self-Induced Oscillations 

The instability discussed in Section 13.21 is induced by the robot itself. Before the system 
gets unstable, it will start to show oscillations. This can be seen for instance in Figured 
When K z = 24.3, the control system’s poles will start to have an imaginary component and 
hence oscillate around the desired value from z = 10m downward. In the field of aerospace, 
such oscillations are a well-known problem, and referred to as Pilot-Induced Oscillations 
(PIO). 
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Figure 4: Root locus plots of a continuous model with a delay of T — 0.15s (as approxi¬ 
mated by a second order Pade transfer function), and aerodynamic drag with \pCr>A = 0.5, 
for z = lrn and z = 10m. The plots at different heights are equal in shape, but for each 
point in the plot K z is different. This is indicated for the (approximations to) the values of 
K z at which the poles get an imaginary component (start of oscillations), and where they 
cross the imaginary axis (become unstable). 


There have been several investigations on the automatic, on board detection of self- 
induced oscillations [321 [ffil Eg S3 E!- The two most important properties of self-induced 
oscillations are that (1) there is a phase shift between the observations and control inputs 
in the order of 90 - 180°, and (2) the oscillations are of a significant magnitude. A typical 
detection method involves a Fast Fourier Transform (FFT) of the observations and control 
inputs in order to detect the occurrence of these properties m u\- In this article, in the 
interest of a computationally efficient and straightforward method, we will investigate the 
use of the covariance between the control input u' z and the relative velocity th and the use 
of d z itself. 

As a first test, landings in simulation will be performed. For the simulations, a ZOH- 
model is assumed with a time step T = 0.03s, \pCr>A = 0.5, a mass of 1 kg, and a delay 
of At = 0.15s (as in Section EZ0 The left part of Figure ED shows a landing for I\ z = 20, 
Zo = 10, v z o = —2, c 1 2 = 0.2, no wind. It shows that the landing goes smoothly, until the 
very end at ~ 0.7m at which point the robot’s height starts to oscillate. 

The right part of Figure [5] shows values that can be observed by the robot itself, zooming 
in at the few seconds before landing. The magenta dashed line is tD z . Toward the end of 
the landing i9 z starts deviating significantly from its reference point —0.2 (indicated with 
the dotted black line). It even obtains positive values, which means that the robot goes 
up. This itself could be used for the detection of instability, were it not that a wind gust 
can also move the robot upward, causing a positive d z at a larger height. The figure also 

1 Please note that while the theoretical analysis involved a linearization of the observation function, the 

simulation employs the nonlinear observation function (and hence can validate some of the conclusions 
from Section 0. 
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Figure 5: Left: Landing in simulation with K z = 20, c 2 = 0.2. The red line is the height, 
the black line is the point at which the onset of oscillations is detected with the robot’s 
observables. Right:. Three of the robot’s observables: d z (magenta, dashed), the thrust 
u' z (blue, solid), and their covariance cov(u' z ,'& z ) (black, dashed-dotted). 

shows the thrust u' z (blue solid line) and the covariance of u' z and d z over the previous 20 
time steps (black dashed-dotted line). During the smooth part of the trajectory cov(u' z , i9 z ) 
is a small negative number: when the robot descends too fast (■$, too negative) the robot 
thrusts more up, and vice versa. 

Towards the end of the landing the observable i9 z starts to change too quickly for the 
control system and co v(u z ,i9 z ) becomes positive: if the robot descends too fast (i9 z too 
negative) the robot will thrust less up, and thus it starts inducing oscillations. Both i} z 
and cov(u{, th) can be used for detecting the onset of oscillations. The black line in the 
left part of Figure 0 shows the first point during the landing at which both variables are 
higher than 0.01. 

In order to test the hypothesis that the onset of oscillations is related to the magni¬ 
tude of the gain and the height of the robot, simulation runs were performed for gains 
K z G {10,30,50} and for various types of wind u win d G {—3, —2.5,..., 2.5,3} m/s. Fur¬ 
thermore, z 0 = 20, v z0 = —1, c 2 = 0.05. The left part of Figure EH shows the results of 
this experiment, with the different K z on the x-axis and the corresponding z at which 
self-induced oscillations are detected on the y- axis. Each K z additionally is represented 
with a different marker, and the markers are colored according to the wind velocity, from 
blue (—3m/s) to red (3m/s). The black dashed line is a linear least-squares fit through the 
points, with as parameters: z = 0.04 K z — 0.1. 

There are three main observations. First, as the linear fit shows, higher gains result 
in oscillations further away from the landing surface. Hence, a fixed gain will result in 
self-induced oscillations around a certain height. The detection of these oscillations can be 
used for instance for triggering landing responses (such as the leg extension by fruitflies). 
Second, a considerable wind difference between —3 to 3 m/s leads to a - for MAVs - 
relatively small difference of ~ 0.50m in the height at which self-induced oscillations occur. 
Third, the colors show that there is a slight correlation between the detection height and the 
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Figure 6: Left: Oscillation detection heights for different gains K z for different wind 
velocities. The color of the marker represents the wind velocity, ranging from —3m/s 
(blue) to 3m/s (red). The black dashed line is a linear fit through the detection heights. 
Right: Results for experiments with a P/-controller, with I~ = 1 as the I- gain. 

wind velocity. Analyzing the landings with wind shows that for the higher wind velocities, 
the control system has more difficulties tracking $*, leading to a steady-state error. This 
suggests that introducing an integrator term could resolve this problem. The results of a 
P/-controller, with K z as the P-gain and P = 1 as the /-gain, are shown in the right part 
of Figure O Since the P/-controllcr can cancel the steady-state error, the reference is 
tracked much better. As expected, self-induced oscillations still occur, and the detection 
heights are much closer together. Moreover, inspection of the marker colors shows that 
there is no obvious relation anymore between the wind velocity and the detection height. 
In order to keep in line with the system that was analyzed theoretically, for the remainder 
of the article we will focus on the P-controller setup with gain K z . 

4.2 Wind Gusts and Actuator Efficiency 

In the previous subsection, the feasibility of detecting self-induced oscillations for distance 
estimation was shown. However, two factors were not modeled that could have a significant 
influence on the practicality of the approach. 

First, during a typical landing outdoors the wind is not constant. The wind can vary 
and sudden wind gusts can occur. In order to study a rather extreme case, wind gusts are 
added to the simulation in the form of a sine function: 

%ust = Wsm(at), (34) 

where W determines the magnitude of the gusts and a the period. 

Second, as mentioned in Subsection 12.21 the mapping from command signal u" to the 
actual force in Newton u' z depends on the actuator effectiveness. In flying robots such as 
rotorcraft or flapping wing vehicles this actuator effectiveness, in turn, depends on the air 
flow. Here, a rotorcraft is assumed and the relation between actuator effectiveness and 
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Figure 7: Left: Landing in gusty wind conditions (ranging from —4 to 4m/s in 7r seconds. 
The red line is the height over time, the black line the detection height. Right: Oscillation 
detection heights for different gains K z for different wind velocities in gusty wind conditions. 
The color of the marker represents the wind velocity, ranging from — 3m/s (blue) to 3m/s 
(red). The black dashed line is a linear fit through the detection heights. 


air flow is modeled according to findings on propellors in |4Tj. Specifically, the following 
formula was used: 

u' z •<— max{)i( — bv air u' z — cu air , 0}, (35) 

so that both the offset and slope of the efficiency change over different air velocities. 

The left part of Figure [7] shows a landing with a rather extreme scenario, in which the 
average wind velocity is Om/s, W — 4, ,o = 1, b — 0.5, and c = 0.5. This means that per 
every 7r seconds, the wind varies from —4 to 4m/s! For this gusty scenario, the threshold 
of cov(u' z ,i!) z ) is set slightly higher, to 0.1. The red line in the figure is the height over 
time, the black line indicates the height at which self-induced oscillations are detected. 
The right part shows the results over many landings with the same settings but different 
wind velocities, from —3 (blue) to 3m/s (red) - to which the gusts are added during each 
landing. The black dashed line is a linear fit through all the detection heights and gains 
K z . 

The results show that despite the very gusty conditions and varying actuator efficiency, 
there is still a positive linear relation between z and K z . The uncertainty is slightly higher 
though than for a constant wind and actuator efficiency, especially for larger gains / heights. 

5 Adaptive Gain Control 

5.1 Distance Estimation in Hover 

In principle, the detection of self-induced oscillations could be sufficient for flying robots 
(or insects) to determine when they are close to their landing target. This allows for the 
triggering of a final landing procedure and hence is behaviorally very relevant. However, 
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Figure 8: Results of simulated hover tests. The robot is placed at different heights with 
K z = 50 with i?* = 0. The gain is then adapted with co v(u' z ,'d z )* = 0.05. x-markers show 
the height z and K z at the first time instant in which co v{u' z , , & z ) > 0.05. The markers are 
color-coded for the wind (blue = —3, red = 3). The black dashed line is a linear fit. 


it would be extremely useful if the robot could determine height at any distance to the 
landing surface. 

Surprisingly, the proposed stability-based strategy for height estimation does not re¬ 
quire the robot to land! Eq. EH] (derived for the vacuum condition) does not depend on 
v z , nor on c 2 . This suggests a strategy for determining the height around hover: the robot 
can set = 0 and change its gain K z so that the control loop starts exhibiting small 
self-induced oscillations. The gain at such a point is then a stand-in for the height. 

Specifically, a control law can regulate co v(u' z , $ z ) by continuously adapting the gain K z . 
This leads to the following adaptive gain control setup. An inner loop uses u z = K z {d* z —d z ), 
while an additional outer loop controls K z : 

K z (t ) = K' z (t) + (P K' z (t))e cov (t), (36) 

K(t) <r- K z (t - T) + (/ K' z )e coy (t - T), (37) 

e co v(t) = co v(u' z , fi z )* - co v(u' z , (38) 

where P, I G [0,1] are a proportional and integral gain for the outer loop control, both 
relative to the current K' z . The reason for this is that larger heights involve higher speeds, 
so that K z should be changed more quickly for larger K' z than for smaller I\' z . 

This strategy has been tested in simulation for z G {5,10,..., 45, 50} and iWd G 
{—3, —2,..., 2, 3}. The simulation stops when |e cov | < 0.005. Figure El shows for each 
simulation the last gain and height with a marker. The markers are color-coded for the 
wind from blue (—3 m/s) to red (3 m/s). The black dashed line is a linear least squares 
fit. 

The main result is that this method indeed gives an approximately linear relationship 
between z and K z , with as linear fit: z = 0.06 K z — 0.2. Although the wind influences the 
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hover altitude (since it is not accounted for in the initial thrust) it again hardly influences 
the results. Finally, despite the theoretical stability analysis, it may be counter-intuitive 
that the relation between K z and z is linear. However, it is instructive to realize that the 
oscillations at 50m are higher than at 5m (in the order of 0.5m versus 0.1m). 

5.2 Distance Estimation During Landing 

Is there a possibility to land and continuously estimate the distance to the surface? Figure 
0 suggests that there is. If the robot descends while keeping co v(u' z , d z ) at a fixed positive 
value, then the gain K z will directly represent the height z during the landing. This 
strategy of “landing on the edge of oscillation” can use the exact same adaptive gain 
control as explained above, but then with c 2 > 0. 



Figure 9: Results of simulated tests landing on the edge of oscillation. The robot is 
placed at different heights with K z = 50. Then, it starts adaptive gain control with 
co <v(u' z , , & z )* = 0.05. First d* = 0, and when for the first time co v(u' z ,i3 z ) > 0.05, the 
relative velocity set point is changed to d* = 0.05. 

Simulations of this strategy have been performed for z G 5,10,..., 50. The implemen¬ 
tation of the strategy starts with a hover maneuver (c 2 = 0) with K z = 50, I = 0.005, and 
P = 0.15. If in hover the condition cov(u , z ,'d z )* = 0.2 is met, the landing is commenced. 
During landing c 2 = 0.05, co v{u' z ,d z )* = 0.05. 

Figure 0 shows the K z versus z during the landing phase. The markers are color-coded 
for the starting height Zo in the landing phase, from blue (50m) to green (5m). After 
switching from c 2 = 0 to c 2 = 0.05, the gain is often too low for causing instability. The 
outer loop has to compensate for this, which causes K z to vary at the start of the landing 
phase. After that, all K z are linearly related to z, with as linear fit: z = 0.1 K z + 0.1. 
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6 Real-world Experiments 

6.1 Experimental Setup 

The experimental setup is shown in Figure [lOj A Parrot AR drone 2.0 is used, running the 
Parrot firmware on board. This means that the state estimation and control of attitude and 
altitude are performed on board, with the help of the drone’s sensors. The sensors consist of 
an IMU (accelerometers, gyros, magnetometer, pressure sensor), and a downward pointing 
camera and sonar. The firmware can be considered a black box, which takes reference 
angles and velocities as commands. If no command is received, the drone cancels lateral 
velocity and keeps its height. 

On the laptop, the in-house developed software package ‘SmartUAV’ is used to com¬ 
municate with the drone [33]. It receives images from the drone at ~ 15 Hz. These images 
are processed with a rather standard vision pipeline. Maximally 40 corners are detected 
with the method in |37], and these corners are tracked to the next image with the Lucas- 
Kanade optical flow algorithm [29]. The method from [9] is used to estimate i? z from the 
optical flow vectors. Second order terms are omitted, as was done in [1] . It assumes a flat 
landing surface and estimates the spatial gradient of the optical flow for determining i9 z . 
This value is low-pass filtered and then used in a control loop v' z = K z {d* — tb.), where 
v' z is the vertical velocity command sent to the Parrot AR drone. Please note that in the 
real-world experiments, the velocity command is saturated at \v' z \ = 0.75m/s. The program 
further allows to adapt K z according to the adaptive gain control strategy explained in 
Section 15.21 Please note that the offboard vision processing setup introduces a significant 
delay between the capturing of the image and the reception of a command signal from the 
offboard laptop, of ~ 1-1.5s. Experiments were performed in an indoor and an outdoor 



Figure 10: Experimental setup. A Parrot AR drone 2.0 sends its images to an offboard 
laptop, which processes the images and sends back control commands. 
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environment, shown in Figure fill 



Figure 11: Left: Indoor test environment. Right: Outdoor test environment. 

Before showing the results of the various tests, it is worthwhile reflecting on the dif¬ 
ference between the real-world experimental setup and the simulation setup. Most impor¬ 
tantly, the main issue at stake here is whether K z contains information on the height z, by 
performing optical flow control on the basis of the relative vertical velocity d z . A major 
difference with the experiments in simulation is that a velocity v' z is controlled instead of 
u' z and that the firmware uses sensors such as sonar to perform this velocity control. This 
may lead to the impression that scaling information is still necessary. However, the value 
of v' z is only used for determining a covariance (whether the control loop wants to go faster 
up or down), and its value in meters per second does not play a role in the estimation. 
Moreover, in the simulations u' z was controlled directly, demonstrating that the use of v' z 
is not necessary to make the approach work. The experiments verify whether the results 
from theory and simulation also apply to a real-world system with significantly different 
actuation dynamics. 

6.2 Indoor Experiments 

First, it is verified that a fixed-gain constant divergence landing results in instability at 
different heights. The result of this test for K z 6 {1,2,3} (green, blue, and red lines) and 
c 2 = 0.1 is shown in Figure [T21 Dotted lines indicate the first point at which co v(u' z , d 2 ) = 
0.1, set to trigger the drone’s automatic landing procedure. The main observation from the 
figure is that also for the real robot, higher gains lead to an instability at larger heights. 

The second finding to verify is the determination of height around hover. For the 
real robot, cov(u z ,'d z )* = 0.05 is employed, with / = 0.075 and P = 2. During the 
hover experiment, the robot has been sent to different heights, subsequently activating the 
adaptive gain control. In order to get an impression of what the adaptive control does, 
let us first zoom in on the relevant variables during a short part of the flight, shown in 
Figure [T3j In this part of the flight, the drone has just ascended to ~ 4.5m (see top left 
plot), while it has a gain of K z fa 2.2 (bottom left). The drone starts to regulate its 
relative velocity to d* z = 0 (top right), and initially there are no self-induced oscillations as 
witnessed by the cov(u{,r? 2 ) ~ 0 (bottom right). The low covariance makes the gain K z 


19 






Figure 12: Three fixed-gain constant divergence landings of a Parrot AR drone, with set 
point d* = —0.1 and gains K z G {1,2,3} (green solid, blue dashed-dotted, and red dashed 
lines, respectively). The dotted lines show the point at which cov(u' z ,'d z ) = 0.1 indicating 
the onset of self-induced oscillations. Given a fixed gain, this instability depends on the 
height. 


increase over time. This gradually causes d z to show ever larger oscillations, which leads 
to higher cov(it', t? z ). Around t = 174s, the covariance reaches its set point, and the gain 
K' z starts stabilizing around 4.6 in order to keep the covariance around 0.05. The question 
now is at what values K' z will stabilize at other heights, and whether there is indeed a 
linear relationship between z and K z . 

The left part of Figure fl4l shows the result of having the drone hover at a given height 
with co v(u' z , , t} z )* = 0.05. The markers are only shown for time instants at which |e cov | < 
0.005. The dashed line is a linear least squares fit. Although more noisy than in simulation, 
the relation between z and K z seems to be roughly linear - as predicted by the theoretical 
derivation. When using the linear fit z — 0.81 K z + 1.1 as a height estimation function, the 
median absolute error of all observations is \z — z\ = 0.28m, while the mean is 0.40m. 

The third finding to verify is the possibility to land on the edge of oscillation. The 
AR drone is sent to a height of 3.5m. Then, the adaptive gain control is started, at first 
with c 2 = 0. As soon as |e cov | < 0.005, the landing phase starts. If K z < 0.5, the drone’s 
automatic landing procedure is triggered. The right part of Figure [14] shows the results for 
all points where |e cov | < 0.005. Different experiments have different markers and colors. 
The dashed line is a linear fit. Although the results are quite noisy, again a roughly linear 
relation appears, with as linear fit z — 0.67A7, — 0.2. The median absolute error of this fit 
is 0.26m, and the mean absolute error is 0.33m. 
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Figure 13: Signals during a part of the indoor hover experiment. Top left: The height 
over time. Top right: Relative velocity over time. The magenta line is the i? z used for 
control, as determined with computer vision, the yellow line is the th as calculated with 
the help of the sonar readings (not used in control). Bottom left: The gain over time. 
The light green solid line is K z , the dark green dashed line is K' z (see Eq. 1361 and 1371) . 
Bottom right: The dashed-dotted line is the covariance of the control input u' z and $ z , 
the dotted line is the set point for the covariance (cov(u' z , th.)*, see Eq. I38]h 
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Figure 14: Indoor results adaptive gain control. Left: Hover experiments at different 
altitudes. The markers represent ( K~,z ) values for which |e cov | < 0.005. Right: Landing 
on the edge of oscillation. The markers are (K z ,z )~pairs during the landing, while K z 
is regulated so that cov(u , z ,$ z ) = 0.05. Markers are only shown at instants at which 
|e CO v| < 0.005. 
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Figure 15: Outdoor results adaptive gain control. Left: Results of hover experiments 
at different altitudes, with co v(u' z ,i} z )* = 0.05. The markers represent ( K z ,z ) values for 
which cov(u' z , i? z ) G [0.005,0.015]. Right: Results of landing on the edge of oscillation. The 
markers are (K z ,z )~pairs during the landing, while K z is regulated so that cov(w',i? 2 ) = 
0.005. 

6.3 Outdoor Experiments 

A few additional outdoor experiments have been performed to confirm that the method 
also works in an outdoor environment with wind and wind gusts. The right part of Figure 
ttH shows the outdoor experimental setup. The ground surface consisted of small stones, 
which already provide ample texture. Some beach toys were added to the scene to also 
ensure visible texture at larger heights. The left part of Figure [15] shows the results from 
the hover experiment. Just as indoors, there is a clear linear relationship between the 
height z as measured by the sonar and the gain K z as determined with the adaptive gain 
control, with as linear fit z = 1.03 K z + 1.4. The median absolute error in height of the 
fit is 0.29m, the mean is 0.35m. The right part of Figure [15] shows the results from two 
landings on the edge of oscillation, again showing a linear relationship. Videos of some of 
the experiments can be seen at@. 

7 Discussion 

7.1 Flying insects 

The presented strategy forms a novel hypothesis on how flying insects such as fruitflies and 
honeybees can estimate distances, viz. by means of the detection of self-induced oscilla¬ 
tions. In [S] a different, thrust-based strategy for distance estimation was introduced, and 
its implications for landing insects discussed. In part, the implications of these strategies 
overlap. They both explain how flying insects can trigger landing responses at a partic¬ 
ular given distance, such as fruitflies that extend their legs just before landing (e.g., 0). 

^https://www.youtube.com/playlist?list=PL_KSX9G0n2P96xb8JRSo0jLFnLBK0Pj5Y 
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However, there is a significant difference in the strategies, which has its consequences for 
the explanation of flying insects’ behavior. A few notable differences in this respect will 
be highlighted below. 

A first difference is in the possible explanation of the hover phase honeybees exhibit 
just before touchdown [13]. The hover phase could in principle be explained as a maneuver 
that is triggered by a distance measurement - for instance by means of the thrust-based 
strategy of [6J. However, the proposed strategy suggests instead that the hover mode is an 
intrinsic property of optical flow control and is part of the distance estimation itself. As 
the honeybee gets closer to the surface, its control gain will start becoming too high. This 
will result in self-induced oscillations and lead to a situation of zero divergence, i.e., hover. 
In this light, it is interesting to remark the characteristic little ‘bump’ in the trajectory in 
Fig. 0 The same bump can be observed in Fig. 4 in Rj. and in Fig. 5 in [39] (in which 
a honeybee is performing a grazing landing). 

A second difference is that the stability-based strategy explicitly allows for the obser¬ 
vation that a visual stimulus itself can trigger the landing responses - even for tethered 
insects ©Eg. The reason for this is that the strategy depends on the control stability 
and hence the system bandwidth. A standard way to measure the system bandwidth is 
to make a Bode plot that maps the system’s frequency response. At too high frequencies, 
the phase shift between the quickly changing observation and the subsequent control input 
will become larger and the magnitude of the control input lower. If the magnitude of the 
control input is large enough, a purely visual stimulus may be detected as a self-induced 
oscillation. 

7.2 Flying robots 

Concerning flying robots, the novel proposed distance estimation strategy provides a novel 
way to perceive distances with a single vision sensor. This may turn out to be essen¬ 
tial for tiny flying robots such as the Robobee [31], for which any type of sensor is a 
significant payload [36] [T9] l2Tj . For slightly larger drones such as 25g ‘pocket drones’ 
(e.g., [□]), the finding is also immediately relevant. For instance, currently the lightest, 
fully autonomously flying robot is the 20g DelFly Explorer |43j, which uses a 4g stereo 
vision system with onboard processing to avoid obstacles in its environment. The pro¬ 
posed stability-based strategy to distance estimation may provide a bigger potential to 
even lighter, single-camera solutions, also for obstacle avoidance. Even for MAVs in the 
order of ~ 1kg or larger, the method may still be useful. It allows to estimate potentially 
large distances without compromising the MAVs’ payload capability. 

However, in order for stability-based distance estimation to live up to its potential, 
some issues need to be further investigated. The experiments have shown the strategy 
to work in an offboard processing scheme with a Parrot AR drone 2.0 with a median 
height estimation error in the order of ~ 0.3m. As a proof of concept this suffices, but 
it is desirable to improve the accuracy of the measurement. This seems quite possible if 
onboard processing is used, so that higher update frequencies can be achieved (e.g., 60Hz 
with onboard processing of the bottom camera images on the AR drone 2.0). 
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8 Conclusion 


In this paper, a stability-based strategy has been proposed with which a robot can estimate 
distances only on the basis of efference copies and optical flow maneuvers. Theoretical 
analysis of linearized models has shown that there is a linear relation between the fixed 
gain and the height at which instability arises. This analysis has been verified in simulation 
(including non-linearities and disturbing factors such as wind gusts) and on a real robotic 
platform. It has been demonstrated that self-induced oscillations can be detected by the 
robot and used to (1) trigger a final landing procedure, (2) determine the height in hover, 
and (3) continuously determine the height during a constant divergence landing. 
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A Appendix: Generalization to other movement di¬ 
rections 

In Section H] it was shown that given a gain K z , the stability of the vertical control loop 
with fi z depends on the height z. In this appendix the case is studied in which the camera 
still looks down in the ^-direction, but the movement is in the horizontal x-direction. 
The control is thus based on the relative velocity fi x = It will first be shown that 
the Ending also generalizes to such horizontal optical flow control (Subsection IA.1I1 . The 
stability-based method then does not require vertical motion. Subsequently, a thrust-based 
method is investigated for horizontal motion (Subsect ion IA.2p . The thrust-based method 
always requires both horizontal and vertical motion. 


A.l Stability-based Method 

Let us start from the observation: 


V = fix = 


(39) 


Since both the horizontal and vertical axes are involved, the state space model will have 
four variables. Linearizing the state space model gives: 


A y(t) = [-J 6 0 0 i] [A z(t) A v z (t) A x(t) A v x (t)\ 
so that the state space model matrices are: 


(40) 
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C = [at 0 0 7 ] ^ = [0]. 


(41) 


Again, the discretized system is studied, which has the following state space model 
matrices corresponding to the continuous ones in Eq. 
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where T is the discrete time step in seconds. The transfer f un ction of the open loop system 
can be determined to be: 


G(w) = C(wl - $)“T, 


(43) 


T 

z(w — 1) ’ 


(44) 


where we use w as the Z-transform variable, since z already represents height. The feedback 
transfer function is: 


G(w) 


K X T 


z(w 


1 ) 


K X T 

z(w—l) 


+ 1 


(45) 


Eq. 1451 shows that, given a gain K x and time step T, the dynamics and stability of the 
system depend on the height z. 

Setting w = — 1 in the denominator and equating it with 0 gives: 


— 2z 



0 


(46) 


, which leads to an equation that expresses the height at which the control system will 
become unstable in terms of the unstable gain value K x : 

* = \k x T, (47) 

which is exactly the same formula as for the vertical movements studied in Section [H 


A.2 Thrust-based Method 


It is instructive to also look at what a possible thrust-based method would look like for 
the horizontal motion direction. For this, let us again start from the observation d x , which 
is differentiated over time: 


• a x v x v z a x n n 

d — — - —- = — — d i9 

u X o U X U Z1 

z Z A z 


(48) 


leading to: 

- ( 49 ) 

{$x + tixtiz) 

Eq. [49] shows that the height can be estimated with the help of a horizontal accelerometer 
measurement and the observables $ x , and the time derivative d x (as was done in [23]). 

If the relative velocity i9 x (sometimes referred to as ‘ventral flow’) is kept constant, then 
d x = 0, and: 

2 ^ or- < 50 > 


where in a vacuum environment, a x can be assumed to be u x . A constant ventral flow and 
divergence landing would then again lead to a way to estimate z, as u x is known and {} x 
and d z can be assumed equal to constants. 
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Interestingly, also this horizontal motion case does not allow for distance estimation 
while staying at the same height ($ 2 = 0). The reason for this can be seen when rearranging 

Eq. EDJ 

a x = z(i9 x tf z ), (51) 

which shows that $ z = 0 implies a x = 0. Indeed, if $ z = 0, 'd x can only be constant if there 
is no horizontal acceleration. 
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